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ABSTRACT 

We explore millimeter line diagnostics of an obscuring molecular torus mod- 
eled by a hydrodynamic simulation with three-dimensional non-local thermody- 
O ! namic equilibrium (nonLTE) radiative transfer calculations. Based on the re- 

sults of high-resolution hydrodynamic simulation of the molecular torus around 
an active galactic nucleus (AGN), we calculate intensities of HCN and HCO + ro- 
tational lines as two representative high density tracers. The three-dimensional 
radiative transfer calculations shed light on a complicated excitation state in the 
inhomogeneous torus, even though a spatially uniform chemical structure is as- 
sumed. We find that similar transition coefficients for HCN and HCO + rotational 
lines lead to a natural concordance of the level population distributions of these 
molecules and the line ratio -Rhcn/hco+ ^ 1 f° r the same molecular abundance 
value over two orders of magnitude. Our results suggest that HCN must be much 
more abundant than HCO + (t/hcn ^ 10 x y HCO +) in order to obtain a high ratio 
(-Rhcn/hco+ ~ 2) observed in some of the nearby galaxies. There is a remark- 
able dispersion in the relation between integrated intensity and column density, 
indicative of possible shortcomings of HCN(1 — 0) and HCO + (l — 0) lines as high 
density tracers. The internal structures of the inhomogeneous molecular torus 
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down to subparsec scale in external galaxies will be revealed by the forthcoming 
Atacama Large Millimeter /submillimeter Array (ALMA). The three-dimensional 
radiative transfer calculations of molecular lines with high-resolution hydrody- 
namic simulation prove to be a powerful tool to provide a physical basis for 
molecular line diagnostics of the central regions of external galaxies. 

Subject headings: galaxies: ISM — galaxies: active — radio lines: galaxies - 
ISM: molecules — radiative transfer 



1. Introduction 

Molecular gas at the center of an active galaxy is a subject of crucial importance for 
observational and theoretical studies of galaxy formation. Various kinds of molecular lines 
such as 12 CO, HCN, H 12 CO + , CS, and so on have been detected in the central regions of a 
number of nearby Seyfert galaxies (Aalto et al. 1995; Curran et al. 2000; Usero et al. 2004; 
Kohno 2005), and luminous and ultra luminous infrared galaxies (LIRGs and ULIRGs : Gao 
& Solomon 2004b; Nakanishi et al. 2005; Papadopoulos 2007; Gracia-Carpio et al. 2007). 
Highly sensitive millimeter telescopes have enabled an increasing number of observations of 
high density tracer lines such as HCN or H 12 CO + , aiming at examination of the density 
structure of the central regions of external galaxies (Curran et al. 2000; Gao & Solomon 
2004a; Gracia-Carpio et al. 2006). Imaging observations of nearby galaxies show complex 
structures in molecular gas in the galactic centers (radius R < 1 kpc), such as compact cores 
and outer rings (NGC 1068 : Planesas et al. 1991; Jackson et al. 1993; Schinnerer et al. 2000; 
Usero et al. 2004, M 51 : Kohno et al. 1996; Scoville et al. 1998; Matsushita et al. 2004 : see 
for reviews of recent results, Kohno 2005; Garcia-Burillo et al. 2006; Kohno, Nakanishi, & 
Imanishi 2007). The similar profiles of CO and HCN lines, of which optical thicknesses are 
suggested to differ, support a small volume filling factor of high density molecular gas (mist 
model; Solomon et al. 1987). Nguyen-Q-Rieu et al. (1992) concluded from the observation of 
15 external galaxies that molecular gas in the central region is composed of numerous clumps 
smaller than the observational beam size of the IRAM 30m telescope (A# H pbw — 26"). 

Zooming into a smaller scale, the unified model of active galactic nuclei (AGN) postu- 
lates a compact, obscuring dusty molecular torus around the central engine (e.g., Antonucci 
& Miller 1985; hereafter we refer to the molecular gas at the central region R < 100 pc of 
an active galaxy as a "molecular torus"). Current observations have failed to resolve the 
small, dense regions within a compact molecular torus of radius < 100 pc. The dynami- 
cal and thermal properties of the putative molecular torus, its structure and origin remain 
ambiguous. 
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High resolution hydrodynamic simulations predict a highly inhomogeneous structure 
inside the molecular torus, which is characterized by clumps and filaments on a subpc to pc 
scale (Wada & Norman 2001, 2002; Wada & Tomisaka 2005). In such an inhomogeneous gas, 
the energy level population should also become complicated, and interpretation of molecular 
line emissions require detailed and careful modeling and calculations. For example, numer- 
ical simulation of 12 CO rotational lines of Wada & Tomisaka (2005) demonstrated a large 
scatter of "X-factor" , a conversion constant between the 12 CO luminosity and molecular hy- 
drogen column density in an inhomogeneous torus. Inhomogeneous density and temperature 
structures are not the unique factors that generate the intricate excitation conditions. There 
are increasing numbers of arguments concerning the effects of peculiar chemical evolution 
by Xray and UV radiation from a central AGN and accompanying nuclear starburst (Kohno 
et al. 2001; Kohno 2005; Imanishi, Nakanishi & Kohno 2006; Imanishi & Nakanishi 2006; 
Imanishi et al. 2004; Aalto et al. 2007; Meijerink, Spaans & Israel 2007; Kohno, Nakanishi, 
& Imanishi 2007; see for a recent review, Garcia-Burillo et al. 2006 and references therein). 
The observed line intensity is the consequence of the complex coupling of various factors, 
such as hydrodynamic and thermal structures, and chemical abundance distribution. Fur- 
thermore, line intensity of high density tracers, such as HCN, may not necessarily represent 
local density (Gracia-Carpio et al. 2006, 2007; Papadopoulos 2007). 

In this article we study the role of high density tracers by the combination of high reso- 
lution hydrodynamic simulation of an inhomogeneous molecular torus and three-dimensional 
radiative transfer calculations without the assumption of local thermodynamic equilibrium 
(nonLTE). We calculate rotational lines of HCN and H 12 CO + (hereafter we omit the isotope 
number 12) and thoroughly examine the excitation status inside the inhomogeneous molecu- 
lar torus. As the first step, we pursue the nonlinear response of nonLTE level population in 
the uniform chemical abundance torus. The chemical evolution in the torus will be discussed 
separately in the subsequent article. 

It should be noted that our calculations can present a direct prediction of the intensity 
distribution of high density tracer lines. Our radiative transfer results are able to bear di- 
rect comparison with observational data. They can provide not only valuable suggestions 
concerning the role of high density tracers (such as HCN and HCO + lines) in current observa- 
tions, but also the theoretical basis for line diagnostics of the thermal and chemical structure 
of a molecular torus. A compact and inhomogeneous molecular torus at the galactic center 
will be an interesting target for future observational instruments, such as the Atacama Large 
Millimeter /submillimeter Array (ALMA). Three-dimensional simulations of a radiation field 
with high resolution hydrodynamic simulations will be an important tool for deriving the 
astrophysical properties from the current and forthcoming high quality observational data. 
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The organization of this article is as follows: In §2 we describe the methods of hy- 
drodynamical and radiative transfer calculations. We display our results of the excitation 
conditions, the intensity and the line ratio distributions, and the relation with the torus 
structure in §3. We briefly discuss the possible effect of chemical evolution on the line ratio 
of HCN and HCO + , and show implications for future observations with ALMA in §4, and 
summarize our article in §5. 

2. Models & Equations 

2.1. Hydrodynamic Model of an AGN Molecular Torus 

We first performed hydrodynamic modeling and simulations of the interstellar medium 
(ISM) around a supermassive black hole following Wada & Tomisaka (2005). An ordinary 
set of hydrodynamic equations is solved with AUSM (advection upstream splitting method: 
Liou & Steffen 1993). We calculated a 64 2 x 32 pc 3 region with 256 2 x 128 uniform Cartesian 
grid points. We employed "model A" in Wada & Tomisaka (2005), in which the molecular 
torus evolves in a steady gravitational potential composed of a supermassive black hole 
(SMBH) of Mbh = 1O 8 M at the center of the galaxy and the galactic dark halo. These 
potential fields are described as $ ex t = — (27/4) 1 / 2 f 2 /(r 2 + a 2 ) 1 ^ 2 for the former and $smbh = 
— G Mbh / (r 2 + b 2 ) 1 / 2 for the latter, respectively. We took the values of the maximum circular 
velocity v c = 100 km s _1 , the core radius a = 10 pc as for the gravitational potential 
of galactic dark halo $ cxt , and the core radius b = 1 pc for the SMBH potential $smbh, 
respectively. The thermal evolution of the torus was calculated with an empirical formula 
of the cooling function of solar metallicity (Z = 1Z Q ), photoelectric heating by strong 
UV background radiation (Go = 10 in Habing unit), and heating by supernova explosions 
within the torus. The parameters that determine the cooling and heating rates, Z and Go, 
are not expected to significantly affect the thermal structure of the turbulent multi-phase 
medium in a torus, since they are determined not only by cooling/heating rates but also by 
dynamical processes, and energy feedback from star formation (Wada & Norman 2001; Wada, 
Meurer & Norman 2002). The torus model assumed a nuclear starburst, which induces an 
average supernova explosion rate of ~ 0.36 yr" 1 . In this model, supernova explosions take 
place randomly in the disk plane. The dynamical evolution of the SNe blast waves in the 
differential rotation in the inhomogeneous torus was explicitly calculated in the simulation, 
as well as thermal evolution due to radiative cooling. 

The structures of density, temperature, and the geometry of the torus are determined by 
the balance of the external gravitational forces (originating in $ ex t and $smbh), self-gravity 
of the gas, turbulent energy dissipation due to radiative cooling, and feedback from super- 
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nova explosions. Wada & Tomisaka (2005) showed that the molecular torus model is globally 
quasi-steady over a timescale of ~ 10 s yr, although local gravitational instability develops 
a significantly inhomogeneous and clumpy internal structure. In our radiative transfer cal- 
culations, we use density, temperature and velocity structures taken from a snapshot of the 
hydrodynamic simulation results. Strong heating by supernova explosions generates a hot 
(T < 10 6 K) and tenuous (n < 1 cm~ 3 ) atmosphere above the dense torus (Wada & Nor- 
man 2001, 2002). In the hot atmosphere, neither HCN nor HCO + molecules are expected. 
Therefore we remove the hot atmosphere from the hydrodynamic simulation results as an 
input to the radiative transfer calculations. The threshold temperature for removing the 
hot atmosphere is assumed to be T th = 1000 K, after comparing with preliminary radiative 
transfer calculations with lower T t h = 400 K. We found little difference in these results, and 
then confidently chose T t h = 1000 K. Thus the temperature range of the input data is 20 K 
<T< 1000 K, and the maximum density reaches as high as nn 2 <2x 10 6 cm -3 . 



2.2. Radiative Transfer Calculations 

We then calculate HCN and HCO + molecular line intensities using the results of hy- 
drodynamic simulations. The input hydrodynamic data is appropriately smoothed from the 
original high resolution results (256 2 x 128) to 64 2 x 32 grid data. The three-dimensional 
nonLTE radiative transfer scheme is adopted by Wada & Tomisaka (2005) for CO rotational 
lines. We calculate nonLTE population distribution simultaneously with ray-tracing along 
randomly sampled rays toward the outer boundary from each cell (for details, see Hogerheijde 
k van der Tak 2000). 

The nonLTE level population n j of the energy level J (where J is the rotational quantum 
number) is calculated by solving the statistical equilibrium rate equation, 

nj^Rjj, = ^{nj'Rj'j), (1) 
jlj^ _ f Ajj> + BjjiJ + Cjji J > J', ^ 

I BjjiJ + Cjji J < J', 

where Ajji, Bjji represent the Einstein's coefficients from the energy level J to J', and Cjji 
is the collisional transition rate per unit time, respectively. In equation (1), the left hand 
side is the outgoing rate from level J under consideration, and the right hand side is the 
incoming transition rate to level J. The average intensity J is calculated from the specific 
intensity I u , 

J = ^- I I i u( j)(v)dvdn, (3) 

4"7T ./ If) 
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where <f>(y) is the normalized absorption coefficient profile (see below). The collisional tran- 
sition rate Cjji is described as 



where nx represents the number density of collision partner, and 7jj/ is the corresponding 
collisional excitation (or de-excitation) coefficient. In this paper we replace nx by the number 
density of molecular hydrogen nn 2 - 

We calculate the specific intensities l v by solving the equation of radiative transfer, 



where S v denotes the source function. Radiatively induced transition rates are calculated 
from J, which is obtained by averaging the specific intensities over the sampling rays coming 
to each grid cell based on the Monte-Carlo method. Rate equations (1), (2), and transfer 
equation (5) are iteratively solved until both the level population and radiation field converge 
self-consistently. The convergence speed is improved by a version of Accelerated Lambda 
Iteration (ALI, Hogerheijde & van der Tak 2000) to be optimizedly tuned for our scalar- 
parallel computer. At the outer boundary we impose the radiation field to be identical to 
the Cosmic Microwave Background radiation (CMB). 

As for transition coefficients (Ajj/, Bjj/, and 7jj'), we adopt the database of Leiden 
University (Schoier et al. 2005). 1 In order to obtain accurate level population distributions, 
we solve the rate equation for < J < J max with the maximum energy level J max = 10 
for both HCN and HCO + lines. The number of sampling rays for J ranges from ~ 300 to 
~ 900 for each grid point to achieve a convergence level \Anj/nj\i < 10 -6 for all J levels 
(Anj\i — n l j — n l J x is defined as the maximum of the difference between rij of the i-th and 
(i — l)-th iterations on the same grid). For the low energy level J (J < 7), the degree of 
convergence progressively improves to Arij/nj < 10~ 10 . 

The model torus has a supersonic turbulent velocity field of a large dispersion Av < 
50 km s _1 , in comparison with the thermal velocity 0.2 km s" 1 < c s < 6.4 km s" 1 . This 
turbulent velocity influences the degree of "overlapping" of the line which accounts for the 
optical thickness. We assume microturbulence in the absorption coefficient profile 4>{y), 




(4) 



x 




(5) 



(f>{v)dv 



Av^/tx 



1 



cxp 



{Avf 



(6) 




1, 



(7) 



LAMBDA : http://www.strw.leidenuniv.nl/~moldata/ 
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where u is the line center frequency measured in the rest frame of each fluid element, 
and Au = (vo/c) ■ or turb is the effective Doppler width described by (T turb , the velocity of 
microturbulence. The velocity structure of the torus is reflected via the profile <f>(v) in 
equation (6). In order to determine a reasonable value for o" tU rb) we survey the parameter 
space ranging between 1 km s _1 < o" tur b < 50 km s" 1 . After confirming the convergence 
stability of the solution over a wide range of parameters 10 km s" 1 < <7 tU rb ^ 50 km s -1 , we 
choose the value of <7 tur b = 20 km s" 1 in the following calculations. 

NonLTE population calculation starts with solution of the rate equation (2) in the 
optically thin limit only with background radiation J = B u (Tqmb) (where B V (T) denotes 
the Planck function) throughout the torus. We examine the effects of the initial guess of the 
population distribution on the final solution by starting with an alternative initial condition. 
As for an opposite extreme alternative, we select the Boltzmann level distribution with 
the local kinetic temperature of each grid for a starting point. We confirm that these two 
initial conditions result in the same solution within errors of the order \Anj/rij\ < 10~ 5 . 
Convergence speed of the radiative transfer calculation is about 2 to 3 times faster for the 
optically thin initial condition than the LTE for both HCN and HCO + lines. We conclude 
that both initial conditions give correct solutions, and then we adopt the optically thin initial 
condition in our radiation calculations because of its faster convergence to the final solution. 

In order to focus on the effects of inhomogeneous torus structure on the excitation 
conditions and emergent line intensities, we assume spatially uniform abundances of emitting 
molecules (HCN and HCO + ). The fractional abundance y = n mo i/nH 2 is assumed to be 
10~ n < y < 10~ 7 . For the same reason, we ignore the radiative pumping by continuum 
emission from AGN and/or nuclear starbursts. The non-uniform abundance distribution 
due to AGN and nuclear starburst will be separately discussed in paper II (Yamada et al., 
in prep). 



3. Results 

3.1. Intensity Distribution 

The left panel of Figure 1 displays the integrated intensity distribution of HCN(1 — 0) 
line of a face-on torus. In this panel we adopt the fiducial value of ?/hcn = 2 x 10~ 9 
for the molecular abundance. The panel exhibits a significantly inhomogeneous intensity 
distribution which reflects the highly inhomogeneous structure inside the torus (see Figs. 
1 and 2 of Wada & Tomisaka 2005 for the density structure). Since level population is 
independent of the viewing angle, hereafter we use the face-on data in the following analysis 
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unless otherwise stated. 

The right panel of Figure 1 displays the line ratio distribution of HCN(1 — 0) and 
HCO+(1-0), R Kcm 

-0)/HCO + (1-0) — 

(jT b (RCN)dv)/(jT b (RCO + )dv) = -Rhcn/hco + (ft is 
the brightness temperature) for the same molecular abundance ?/hcn = 2/hco+ = 2 x 1CT 9 . 
The panel shows that in spite of the significant inhomogeneity of the intensity distributions, 
the line ratio is restricted in a narrow range from ~ 0.2 to ~ 1.2 (the probability distribution 
function of -Rhcn/hco+ nas a snar P peak around the median ~ 0.6). The similarity of 
HCN(1 — 0) and HCO + (l — 0) intensities arises from the almost identical values of rotational 
constant B and permanent electric dipole moment \x e for these molecules, which determine 
the energy levels and Einstein's A and B coefficients for pure rotational transitions (see 
Table 1). The coefficients in Table 1 are calculated by the equations 

Ej = BhJ(J + l), (8) 



h 



A, hJ . x = 16- B 6 J 6 Bj,j_ 1: (10) 



where Ej is the energy measured from the ground state, h is the Planck constant, and c 
is the velocity of light. Besides these coefficients, the critical density for LTE distribution, 
n crit = Aj > ji/ r yj i ji is different between HCN and HCO + molecules. However, the amount of 
gas mass of densities larger than the critical density of HCN(1 — 0) (n crit ~ 10 6 cm" 3 ) is 
~ 3.6 x 1O 5 M , and that larger than critical density of HCO + (l - 0) (n crit ~ 10 5 cm -3 ) is 
~ 1.3 x 1O 6 M , which agree within a factor of order unity. In the optically thin limit, line 
ratio straightforwardly traces the mass ratio, which becomes about unity in our simulation. 
In our calculations average optical thickness over the entire field of view is at most (r ) = 2 
(see §3.3 below), and thus dependence of line ratio on critical densities is weak. 

The overall tendency -Rhcn/hco+ — ^(-0 does not depend on the molecular abundance 
y. In Figure 2 we present the integrated intensities averaged over the field of view (64 x 
64 pc 2 ) as a function of y for both HCN(1 — 0) and HCO + (l — 0) lines. Figure 2 shows 
that (i)HCN(i-o)(y) (-Ohco+(i-o)(z/) f° r over ^ orders of magnitude of y. This inequality 
(-Ohcn(i-o)(j/)/ ( j )hco+(i-o)(i/) < 1 implies that in order to obtain a high ratio i? H CN/Hco+ ~ 
2 observed in some nearby galaxies (e.g., NGC 1068, Usero et al. 2004), HCN molecules 
should be much more abundant than HCO + molecules (j/hcn 2/hco+)- For example, if 
we take 2/ H co+ = 2 x 10 -8 , Figure 2 means that yncN should be as large as > 10~ 7 (or 
2/hcn > 10 x y H co+) for a hi g h ratio ^hcn/hco+ - 2 - 
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3.2. The Line Ratio -Rhcn/hco+ i n the Bright Regions 

The left panel of Figure 1 shows that the torus has many bright spots of several parsec 
in size in the integrated intensity distribution. The line ratio distribution has a smoother 
structure with -Rhcn/hco+ ~ 1 (the right panel of Fig. 1). Unfortunately the observational 
beam size of present instruments is too large to resolve such small structures inside a com- 
pact torus (for instance, A# H pbw = 2" ~ 6" for NMA and Rainbow interferometers typically 
correspond to several hundred pc except for Local Group members). Observational estima- 
tion of line ratio would be practically governed by the ratio of the bright regions within the 
unresolved torus. In the following paragraph, we examine the line ratio -Rhcn/hco+ hi the 
bright region, which is expected to be close to the observational estimation in the unresolved 
molecular gas. 

We show two-dimensional probability distribution function of the integrated intensity 
I(x,y) and the ratio -Rhcn/hco+^j v) = ^hcn(z, v)/Ikco+( x , v) in Figure 3 (where I(x,y) is 
the integrated intensity at the position (x,y) in the field of view). Three panels in Figure 3 
are for different molecular abundances, y = 2 x 10~ 10 , 2 x 10~ 9 , and 2 x 10 -8 , respectively. 
Figure 3 shows the line ratio -Rhcn/hco+ ~ 1 f° r almost all the pixels irrespective of the 
integrated intensity I(x,y). When y is very small (= 2 x 1CT 10 ), pixels of high ratio up to 
-Rhcn/hco+ < 1.4 remain at the bright end. As y approaches a realistic abundance value 
y = 2 x 1CT 8 inferred by observational results (e.g. NGC 1068, Usero et al. 2004), however, 
the line ratio at the bright end falls below -Rhcn/hco+ ~ 0-8- The lower -Rhcn/hco+ towards 
the brighter points in the case of y = 2 x 10~ 8 may indicate the possible underestimation 
of the line ratio -Rhcn/hco+ i R the current observations of spatially unresolved tori. This 
trend in Fig. 3 might even strengthen the necessity for abundant HCN over HCO + to obtain 
a high ratio (-Rhcn/hco+ ~ 2). On the other hand, it might be possible that different spatial 
distributions of HCN and HCO + molecules loosen this requirement, but this topic is out of 
this paper and will be discussed in a subsequent paper. 



mol. 


2 J B[GHz] 


He 


A w [Hz] 


Bio 


HCN 
HCO+ 


88.63 
89.19 


2.99 xlO- 18 
3.93 xlO- 18 


2.40 xlO" 5 
4.25 xlO" 5 


2.35 x 10 9 
4.25 x 10 9 



Table 1: Fundamental parameters of HCN and HCO + pure rotational transitions. The value 
2B is twice the rotation constant and is identical to v w . Electric dipole moment \i e is 
measured in cgs units, A w measures the number of spontaneous transitions per unit time, 
Bio is the stimulated radiative transition coefficient per unit time and per unit intensity, 
respectively. 



3.3. The Population Distribution in the Radiative Transfer Calculations 

In this subsection we demonstrate the population distribution of our nonLTE calculation 
results (the basic characteristics of multi-level population distribution in statistical equilib- 
rium in the optically thin limit are summarized in Appendix). In Figure 4, we plot the ratio 
ni/n against density n Ha for three models with y = 2 x 1CT 10 , 2 x 1CT 9 , and 2 x 1CT 8 . In 
Figure 4, the level population in the optically thin limit is also shown with lines. Figure 4 
exhibits two distinct features of the ratios ni/rio of radiative transfer calculations compared 
with optically thin cases: 1) larger values of (ni/no) > (^i/^o)thin on the low density side 
( n n 2 ~ 104 cm" 3 ), and 2) smaller peak value of (n 1 /n ) pcak < (™i/™o) P cak,thin at n H , 2 ~ n crit . 
The departure from optically thin distribution escalates with the molecular abundance y. 
These features on the diagram can be elucidated by the increase in the radiation-induced 
transition rate (oc J) by the propagation of line emission in the torus. 

The effect of frequent radiation-induced transitions appears on the diagrams as 1) the 
broad dispersion in the n\jn§ distribution on the low density side (10 2 cm -3 < nn 2 J$ 10 4 
cm -3 ), and 2) the decrease in the local peak value of ui/uq at the density n ~ n cr it. In the 
low density regime (n <C n cr - lt ), collisional transitions are negligible and the rate equation 
becomes 

n B 01 J = ni^io + ni-B 10 J, (11) 
and the ratio Ui/uq is written as 

rh _ gi/goB w J 
n A w + B w J 

Equation (12) means that ni/n increases with the average intensity J, and can become as 
large as n\/no < gi/go = 3. Therefore ni/no can take a value ranging from 0.67 (when 
J is equal to the background radiation, Jcmb; see Eq. [Al]) to 3 according to J. In 
the denser regime, the decrement of the local peak value of rii/rio at nn 2 ~ n cr[t is also 
explained by the similar argument in Appendix A.l as follows: A large radiative excitation 
rate (oc J±o = 1/(Att) f f I u (j){yio)dudVL) increases the ratio rii/rio, and then the excitation 
rate from the level J = 1 (not from the ground level J = 0) increases as well. The increase 
in the excitation rate from level J = 1 obstructs the downward cascade J — 2 — > 1, and thus 
leads to the reduction of particle accumulation at J = 1 (see Appendix A.l for details). 

Increase in molecular abundance y raises the local emission rate (j v = a v S v in Eq. [5] 
increases as oc n mol n H2 oc y aio \ n \^)-, and then the average J increases with y as well (Eq. [3] 
and [5]). Hence the radiation- induced transition rate due to J increases with y. The increase 
in J io reduces the peak value of n\jno at nn 2 ~ n cr i t . The dispersion of ni/n in the low 
density regime and the degree of reduction in the peak value of ni/n at n Ha ~ n crit become 



(12) 
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larger for a larger y in Figure 4. This behavior also supports the importance of radiation 
induced transitions (oc J) in determining the excitation conditions in the inhomogeneous 
torus. 

In addition to the increased radiation-induced transitions, the local photon trapping can 
also affect the distribution of ni/n for a large y case. The effect of the local photon trapping 
can be described in terms of the escape probability (3 (Goldreich & Kwan 1974; Peraiah 2002). 
The local photon trapping reduces the effective critical density as n crit (eff.) « f3A U) / r ) V) . The 
escape probability (3 is described as f3 — (1 — e -3r °)/3To for a gas slab and f3 — (1 — e -T(, )/T 
for a spherical gas (e.g, Peraiah 2002). Though our model torus is neither a slab nor a 
sphere, both formula provide limitting values of (3 ~ 1 for small r (r < 1) and (3 ~ r -1 
for a large r (r 1). The average optical thickness in the field of view is estimated to 
be (r ) = 2.00 for the face-on data at the line center of HCN in our simulations for a large 
y — 2 x 10 -8 result. The average optical thickness of order of unity in our calculations results 
in the similar value of (3, so that only a negligible effect due to (3 is expected on the diagrams 
of Figure 4. 

3.4. The Effects of Inhomogeneous Structure on Optical Thickness and 

Intensity Distributions 

In the previous subsections we have examined the statistical characteristics of the ex- 
citation conditions in the torus. In order to examine the connection between the resultant 
emission and the input torus properties, in Figure 5 we present the integrated intensity as 
a function of column density Nn 2 . Figure 5 shows that while the integrated intensity / is 
tightly correlated with Nn 2 in the low 7Vh 2 regime (A^h 2 < 10 23 cm -2 ), the dispersion of / 
appears at the larger N^ 2 regime (iV"H 2 > 10 23 cm -2 ). The origin of this dispersion is further 
investigated in the relations of integrated intensity, optical thickness, and column density 
(Fig. 6). Figure 6 (a) shows a dispersion in the optical thickness t = f a u ds at large 
column density (N H2 > 10 23 cm -2 ). This dispersion is obviously produced by the diversity 
of the absorption coefficients a u (rin 2 ,T) due to the nonLTE level population in the inhomo- 
geneous torus. The drastic increase in integrated intensity appearing at r < (Fig. 6 [b]) 
owes the stimulated emission by the population inversion (Appendix §A.2). 2 In the two 



2 In our simulation, the gain factor e -ro due to population inversion is at most a several times of 10 (see 
Fig. 6 [6]). This gain factor might be underestimated by the adopted microturbulence, oturb = 20 km s -1 in 
counting for to in a grid (Eq. 6). Higher resolution simulation might find even stronger masar spots amplified 
by the population inversion and velocity coherence, though this kind of calculation is computationally too 
hard a task at present. 
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panels of Figure 6, crosses indicate the values of the average integrated intensity ((/) = 42.0 
K km sec -1 ), the average optical thickness ((t ) = 0.37), and the average column density 
((-/Vh 2 ) = 9.27 x 10 22 cm -2 ) in the face-on field of view. Our three-dimensional nonLTE 
radiative transfer calculations reveal the large dispersions around the average values of /, 
To, and Nn 2 . These results imply that high angular resolution observations of ALMA are of 
crucial importance to study the structure of the torus. 

The spatial correlations of optical thickness and integrated intensity also show a wide 
variety in our results. In Figure 7, we plot the spatial distributions of the integrated intensity 
of HCN(1 — 0) and the optical thickness To. The disaccordance of the integrated intensity 
distribution (color map) and the optically thick region (to > 1) is obvious in this panel. 
This disaccordance is reasonably well explained by the different dependence of r and I on 
density and the line-of-sight structure of the molecular torus as follows. 

In three panels of Figure 8, the distributions of Ar (the optical thickness per single 
grid) and the density along three lines of sight are displayed. The locations of three lines of 
sight are indicated in Figure 7. The three lines of sight represent different circumstances in 
the inhomogeneous torus. In panel (a) emission from a clump of % 2 < 10 3 cm" 3 is visible 
through a tenuous and optically thin ambient medium, in panel (6) a tenuous atmosphere 
encompasses a large scale height in the outer torus, and in panel (c) the stimulated emission 
due to the population inversion inside a high density region around z = plane dominates 
the intensity. Since the local emission rate j v = a v S v (oc nHCN^H 2 n H 2 ) an< ^ the optical 
thickness r (oc Nn 2 oc nn 2 ) are differently dependent on density, the integrated intensity can 
be strong even if To is not large in an inhomogeneous torus (Fig. 8 a and c). When population 
inversion {a v < 0) occurs (in panel c), the integrated intensity can be much larger than that 
simply estimated from the column density (Fig. 5). This is because while the intensity rises 
rapidly in proportion to el Ar °l, the optical thickness r = J a ds only weakly increases due 
to the negative a in the integrand. 

Our results show that even if the chemical abundance distribution is uniform, the dis- 
tribution of the intensity and the optical thickness can significantly differ. Hence it would 
be possible that the line ratio reflects neither the ratio of the corresponding amount of the 
gas of density higher than critical density, nor the ratio of the area of the r = 1 surface. 
Conventional arguments about the amount of dense gas or the r = 1 surface area are ap- 
plicable only if the optical thickness is assured to be quite small (tq <C 1) or large (tq ^> 1). 
However, if the ISM is subthermally populated and/or has r ~ 1 (which would be the case 
in the molecular torus), one should take account of the large dispersion of a u , especially of 
population inversion. 



4. Discussion 



We performed three-dimensional radiative transfer calculations with the assumption 
of spatially uniform chemical abundance, and closely examined how the inhomogeneous 
torus structure affects the excitation conditions and the line ratio -Rhcn/hco+- ^ ur results 
in §3 are briefly summarized as follows: 1) the expected line ratio of the same molecular 
abundance y is -Rhcn/hco+ ~ 1 over a wide range of y, and 2) the intensity can be strongly 
affected by the stimulated emission by the population inversion. The arguments in §3.2 
and Figure 3 indicate the smaller line ratio -Rhcn/hco+ ^ 1 i R bright regions, and thus 
the observational estimation of line ratio of the unresolved torus could not be the "true" 
average ratio -Rhcn/hco+ including the faint regions. Even if we consider the inhomogeneous 
thermal structure in the torus, a ratio of -Rhcn/hco+ > 1 is quite difficult to achieve with 
the assumption of spatially uniform chemical abundance (§3.2). 



4.1. HCN and HCO + Line Ratio as a Probe of Chemistry of Molecular Gas 

Molecular line intensities could reflect the chemical evolution of molecular gas under the 
irradiation of Xray and UV emissions from the galactic center. This kind of study began 
with the pioneering work of Kohno (2005) on the "HCN diagram" of nearby Seyfert and 
starburst galaxies. His "HCN diagram" shows a trend of stronger HCN line emission in 
nearby galaxies with little signature of starbursts (-Rhcn/hco+ ~ 2.0 and -Rhcn/co ^ 0.4), 
compared with those with starburst regions (see also Imanishi, Nakanishi & Kohno 2006; 
Imanishi & Nakanishi 2006; Gracia-Carpio et al. 2006). In addition to local Seyfert galaxies, 
a number of other types of galaxies show high ratios -R H cn/hco+ > 1-0 (f° r example, NGC 
4418; Imanishi et al. 2004, UGC 5101 and Mrk 273; Imanishi & Nakanishi 2006; IC 342, 
Maffei 2, NGC 2903, Nguyen-Q-Rieu et al. 1992 ; NGC 1097, NGC 1068, NGC 5194, Kohno 
2005). In this subsection we briefly discuss our radiative transfer results and chemistry. 

Our results show the ratio -Rhcn/hco+ ~ 1 f° r the same molecular abundance y over 2 
orders of magnitude (the right panel of Fig. 1 and Fig. 3). The arguments in §A.l derive the 
reason for the convergence of the ratio ni/riQ of HCN and HCO + in both the extremities of the 
low and high densities in the optically thin limit (Fig. 11). Since the intensity I\q increases as 
yrii, the intensity ratio -Rhcn/hco+ rou ghly agrees with the fractional level population ratio 
rii/n of HCN and HCO + . Hence for the same value of y, line ratio -R H cn/hco+ becomes order 
of unity (Fig. 3). 3 This result means that, in order to obtain a high ratio -Rhcn/hco+ > 1> 



3 The brightness temperature Xj, is proportional to v 2 , and therefore -Rhcn/hco+ * s accordingly reduced 
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HCN should be much more abundant than HCO + (t/hcn 3> Unco+)- Recently Gracia-Carpio 
et al. (2007) found that the high line ratio requires about 10 times larger HCN fractional 
abundance than HCO + from their Large Velocity Gradient (LVG) analysis of J = 3 — 2 and 
J = 1 — transitions of LIRGs and ULIRGs. Our results are approximately consistent with 
theirs. 

Several groups discussed the chemical evolution scenario for a high ratio -Rhcn/hco+- 
Meijerink, Spaans & Israel (2007) recently examined chemical evolution of Xray- and UV- 
irradiated molecular gas with a wide variety of input parameters. They calculated thermal 
and chemical evolutions of plane parallel gas slabs. Their results show that in a photodisso- 
ciation region (PDR) the ratio tends to be high (-Rhcn/hco+ > -0 * n an Xray dominated 
region (XDR) the tendency is reverted (-Rhcn/hco+ < -0> exce Pt for a small number of 
models of high density and strong Xray emission. Though the majority of these kinds of cal- 
culations assumes a simple geometry for molecular gas (e.g., Maloney, Hollenbach & Tielens 
1996; Meijerink & Spaans 2005; Meijerink, Spaans k Israel 2007), PDR and XDR chemistry 
currently does not seem to explain the "HCN diagram" classification (Kohno 2005) of AGN 
and starbust galaxies. As for alternative scenarios, the intense ionization flux from nearby 
supernova remnant shock waves (and accordingly high ionization degree, Lepp & Dalgarno 
1996) and shock chemistry have been proposed (e.g., Nguyen-Q-Rieu et al. 1992), but it 
has been still uncertain whether these ISM chemistry models can find the evolutinonary 
path for yielding 2/hcn/2/hco+ > 1- Besides the chemical abundance models, Aalto et al. 
(2007) discussed the possible role of MIR photons from the dusty torus in determining the 
excitation conditions of the emitting molecules for LIRG/ULIRGs. 

We expect that the molecular torus should have an inhomogeneous structure, and then 
the shielding of high energy (UV and Xray) photons that affects the chemical and thermal 
structures would be significantly different from the those predicted by one-zone or one- 
dimensional cloud models. Three-dimensional chemical abundance effects, such as simple 
models of PDR and XDR will be investigated with the three-dimensional scheme in the 
subsequent paper. 



4.2. Implications for Future Observations 

We mainly examine molecular line emission in millimeter band on the basis of the 
nonLTE radiation transfer simulation so far. If we consider a nearby galaxy at the distance 
D = 20 Mpc, the size of the bright spots in our results (Fig.l), of which radius 1 pc, 



for the same value of ni/n . 



- 15 - 



corresponds to the angular size A9 m 0.03". Therefore internal structures of a torus will 
be resolvable with ALMA in the line intensity distributions. The dispersions of Nr 2 , a u , 
and the intensities presented in Figures 7 and 8 will be revealed in detail by ALMA's high 
angular resolution observation as well. We briefly discuss the properties relevant to future 
observations. 

Equation (6) means that our radiative transfer calculations take into account the turbu- 
lent velocity field inside the molecular torus. Figure 9 demonstrates the line profiles with 10 
pc spacing. We assume a nearby galaxy at the distance D ~ 20 Mpc and take the binning 
size of the velocity component along the line of sight to 16 km s _1 . The line profiles in 
Figure 9 show an extremely complicated structure, reflecting the inhomogeneous structure 
and the turbulent velocity field in the torus. Since a line of sight can pass through more than 
one dense clump, multi peaks appear in the line profile, even if the molecular abundance is 
uniform. Figure 10 presents the average profile over the entire field of view. The compiled 
profile in Figure 10 shows a Gaussian-like structure of single components: this means that 
even if current observational results present Gaussian-like profiles, future high-resolution ob- 
servations will be able to reveal the internal substructures of the compact molecular torus 
(i? < 100 pc). 



5. Summary & Conclusion 

We performed three-dimensional nonLTE radiative transfer calculations for HCN and 
HCO + rotational lines based on high-resolution hydrodynamic simulation of an AGN molec- 
ular torus. An AGN molecular torus is expected to have exceedingly inhomogeneous density 
and temperature structures, and accordingly complicated energy level populations. In this 
article we examined the nonLTE level population distribution in the inhomogeneous molec- 
ular torus and its effects on the line ratio -Rhcn/hco+ w hh the assumption of the spatially 
uniform chemical abundance distribution. The results are summarized as follows: 1) The ra- 
tio of HCN and HCO + rotational lines becomes -R H cn/hco+ ~ 1 f° r a w ^ e ran g e of molecular 
abundance y (10~ n < j/hcn = 2/hco+ — 10 7 )> thus obtaining a high ratio -Rhcn/hco+ ~ 2 
observed in some galaxies (such as NGC 5194 and NGC 1068, Kohno 2005) requires more 
abundant HCN than HCO + ((?/hcn) ^ 10 x (?/hco + ))- 2) The spatial distribution of line ratio 
i? HCN y HCO + is inhomogeneous, reflecting the spatially inhomogeneous structure of the molec- 
ular torus. 3) Inhomogeneity in the structure of the molecular torus generates a dispersion 
around the linear relations between intensity and column density. The stimulated emission 
by population inversion can dominate the integrated intensity where riu 2 ~ ?Vit (panel c of 
Fig. 8). When the gas is subthermally populated and marginally thick (r ~ 1), the line ratio 
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may indicate neither the fraction of high density molecular gas nor the ratio of surface area 
of To = 1 region. Furthermore, our three-dimensional nonLTE calculations demonstrate the 
complex line profile distribution in the synthetic image of the torus (Fig. 9), even with the 
uniform abundance assumption. The forthcoming ALMA's high-resolution observations will 
reveal such complex structures. The compilation of high-resolution hydrodynamic simula- 
tions and three-dimensional radiative transfer calculations will open a powerful way to study 
the compact molecular gas in central regions of external galaxies. 

We are grateful to K. Kohno and K. Nakanishi for their helpful comments on the mil- 
limeter observations, and to M. Imanishi for his comments on IR observations of LIRGs 
and ULIRGs. MY thanks K. Sakamoto and S. Takakuwa for their advice on the analysis 
of radiative transfer results. This research was supported in part by Grant-in-Aid by the 
Ministry of Education, Science, and Culture of Japan (16204012, 17340059, 18026008). 

A. Excitation Temperature Distribution 

In the Appendix we describe the fundamental mechanisms that determine nonLTE level 
population or excitation temperature distribution in the optically thin limit. We especially 
investigate the superthermal level population appearing at n ~ n crit , because the stimulated 
emission due to the superthermal population plays an important role for the line intensities 
of high density tracers. In order to provide a simple description of the populating processes, 
we here focus on the low energy level excitation in terms of the ratio ni/n , compared with 
the two level system. 

A.l. Population Distribution in the Optically Thin Limit 

We first study how the ratio ni/rio behaves according to the increase in the number 
density of collisional excitation partner nn 2 in the optically thin limit. In Figure 11 the 
fractional population ni/n is plotted as a function of density for various temperatures. We 
solved equations (1) and (2) ignoring the average line intensity J except for the cosmic 
microwave background radiation (J = B v (Tqmb) = Jcmb)- In Figure 11, HCO + and HCN 
populations are shown in black and blue lines, respectively. It is obvious in Figure 11 that 
the fractional population ni/n of HCN coincides with that of HCO + at the extremities 
of low (n H2 < 10 2 cm" 3 ) and high (n H2 > 10 7 cm -3 ) densities for the same temperature. 
This fact can be easily understood in terms of the multi-level population of pure rotational 
transitions as follows: 
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In the low density limit, collisional excitation by molecular hydrogen is negligible, and 
level population distribution can be calculated by the balance of radiative excitation by CMB 
photons and spontaneous emission decay (Eq. [11]). Inserting J = Jqmb m equation (11), 
the ratio ni/n becomes constant for both of HCN and HCO + , 

ni _ gi/gogioJcMB _ 67 ( A1 ) 
no A w + Bi J C MB 

In the second equality in equation (Al) we use the fact that Einstein's A and B coefficients 
for HCN and HCO + rotational transitions are almost the same (Table 1). 

On the other hand, in the high density limit, the level population converges to the 
Boltzmann distribution, 

nj gj ( Ej - Ej, 
exp 



nj' g,j> V k B T 



2J+1 



Bh{J{J + 1) - J'(J'+ 1)} 



From Table 1, rotational constants £?(HCN) ~ Z?(HCO + ), and then the Boltzmann distri- 
butions of two lines are in perfect agreement, 

Vlte" 6XP V k^Tj' 

Equations (Al) and (A2) explain the natural concordance of ni/n of HCN(1 — 0) and 
HCO + (l — 0) in both low and high density extremities. 

Figure 11 exhibits a peak of (ni/no) higher than that expected for LTE (71i/tio)lte a ^ 
the intermediate density (10 2 cm -3 < % 2 < 10 7 cm -3 ) for each temperature. This peak 
is a result of a combination of two effects, namely the energy level cascade in a multi-level 
system, and the strong J dependence of critical densities for pure rotational transitions. 

The energy level cascade originates in the selection rule of radiative transitions. Pure 
rotational transitions have a selection rule of A J = ±1. The selection rule significantly 
simplifies the rate equation in the low density and optically thin limit (Yamada et al. 2007), 

njAjj-x = ^2 C o,J' n o- (A3) 
j'>j 

Equation (A3) means that a collisionally excited particle from the ground state should cas- 
cade down all the steps of the energy level ladder with A J = — 1 accompanied by spontaneous 
emission photons. 
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As density increases, collisional de-excitation begins to modify the simple form of rate 
equation (A3). As shown in equations (9) and (10), the Einstein's A coefficient for pure 
rotational transition is strongly dependent on J as Aj j_i oc J 3 . On the other hand, the 
collisional transition constant only weakly depends on J (Goldreich & Kwan 1974; McKee et 
al. 1982). Thus the critical density for LTE n CT \ t (J, J — 1) = Aj 5 j_i/7j 5 j_i sharply increases 
with J (oc J 3 ). Figure 12 displays the fractional level population distributions fj (fj = 
n,j/n mo \, J2j f.J = 1) as a function of J for various densities (n H2 = 10 5 , 10 7 , and 10 9 cm -3 ). 
The overall distribution of fj approaches the Boltzmann distribution (indicated with a solid 
line in Fig. 12) as the density increases. In Figure 12, it is also shown that while fj is close 
to the Boltzmann distribution for low J, fj deviates from the Boltzmann distribution for 
high J. This is another expression for the reason the peak of ni/n$ appears (Fig. 11): for a 
fixed density n, the low J levels that satisfy n ^> n cr i t (J, J — 1) approaches the Boltzmann 
distribution (the solid line in Fig. 12), and the further downward level cascade from the 
corresponding level comes to a halt. However, at higher J levels where the LTE condition 
is not satisfied, the downward cascade by spontaneous emission accumulates the cascading 
particles at the density n ~ n crit . Accumulation of the extra cascaded particles in addition 
to the Boltzmann distribution generates a peak of rii/n shown in Figure 11. 



The source function S v and the absorption coefficient a v from the level J to J — 1 denote 



These equations imply negative S v and a u if 5 pop = njgj_i/rij_igj exceeds unity (popula- 
tion inversion: Rybicki & Lightman 1979). In Figure 13 we plot S pop for HCN(1 — 0) and 
HCO + (l — 0) in the optically thin limit as a function of density and temperature. In our 
torus model, density and temperature range riu 2 <2x 10 6 cm -3 and 20 K < T kin < 1000 
K, respectively. Figure 13 suggests <5 pop is likely to become larger than unity in this density 
and temperature regime. 

When population inversion occurs in the i—th grid (S pop > 1), optical thickness in the 
grid At* = a % u As becomes negative due to negative absorption coefficient a l u . It is apparent 
in terms of the formal solution of radiative transfer I v {t v ) = I u (ti)e~ Tv + S u (l — e~ Tv ) that 
negative optical thickness t u < strongly enhances the specific intensity I u . Since population 
inversion can significantly affect the intensity, it should be treated correctly in line transfer 



A. 2. Population Inversion 



S, 




(A4) 



(A5) 
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calculations. As described in §A.l, one of the origins of population inversion between J 
and J — 1 (or peak of Ui/uq in Fig. 11) is the downward energy level cascade from higher 
levels j (> J) to J. Therefore the maximum energy level J max in the radiative transfer 
calculation has to be as large as possible to obtain the precise level population, especially 
when population inversion is expected. Preliminary radiative transfer calculations reveal the 
lower limit of the maximum energy level should be J max > 8 in our model torus. In this 
paper, we take the upper limit of J max =10. 

Intense population inversion (S pop > 1) occasionally takes place during the iterative 
calculations before the solution of level population meets the final convergence. This sort of 
temporal population inversion sometimes induces numerical divergence of the intensity due 
to a factor e~ Tv . We avoid the temporal numerical divergence using a way similar to that 
of Hogerheijde & van der Tak (2000): we set the numerical upper bounds on J and \t v \. 
Numerical experiments with several sets of upper bounds present an excellent convergence, 
thus confirming the stability of the obtained solutions against the variation of the numerical 
limiters. 
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Fig. 1. — Left panel: integrated flux distribution of HCN(1 — 0) line with yncN = 2 x 1CT 9 
for face-on data. Right panel: the distribution of the line ratio -Rhcn/hco+ with the same 
abundance value of y. In the right panel, the maximum number density of hydrogen molecules 
on all the lines of sight are overplotted with green contours. 



1200 




10" u 10"'* 10* 9 10 s 10" 7 

y 



Fig. 2. — Spatially averaged integrated intensity (I) over the entire field of view as a function 
of molecular abundance y. A line with squares displays (I) of HCO + (l — 0) emission, and 
that with crosses presents (/) of HCN(1 — 0) emission, respectively. 
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Fig. 3. — Combined probability distribution function of the integrated intensity and the 
ratio -Rhcn/hco+- Abscissa denotes the integrated intensity of HCN (1 — 0) line and ordinate 
denotes the ratio of Jhcn(i-o)/-^hco+(i-o) P er pixel, respectively. Three panels correspond to 
different abundances y, 2 x 10~ 10 , 2 x 10 -9 , and 2 x 10~ 8 from left to right. For a small 
value of y, the line ratio -Rhcn/hco+ becomes ~ 1 in all pixels. As y increases, though, the 
number of pixels of low ratio -Rhcn/hco+ — 1 accordingly increases in bright pixels. 
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Fig. 4. — The relative population ni/n Q as a function of density rin 2 for HCN (the left panel) 
and HCO + (the right panel). Lines are those for the optically thin limit (the same as Figure 
13) of which kinetic temperatures are 10, 30, 100, 200, and 1000 K, respectively. The results 
of our radiative transfer calculations for three values of y are presented with dots of different 
colors: black dots are of y = 2 X 10~ 8 , pink dots are of y = 2 X 10 -9 , and dark red dots are 
of y = 2 x 10~ 10 . 
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Fig. 5. — The integrated intensity as a function of the column density iV H2 of the face-on 
data. When Nu 2 is larger than ~ 10 23 cm" 2 , the dispersion of the integrated intensity / 
becomes accordingly larger. A cross at A H2 = 9.27 x 10 22 cm -2 and / = 42.0 K km s -1 
indicates the average column density and the average integrated intensity in the field of view. 




Fig. 6. — The distribution of optical thickness r and column density A H2 (panel a), and of 
optical thickness and integrated intensity (panel b). The data are taken from the results of 
HCN(l-O) line calculation with y = 2xl0" 9 . Crosses at (iV H2 ,r ) = (9.27xl0 22 [cm" 2 ], 0.37) 
in panel (a) and (t ,J) = (0.37,42.0 [K km s" 1 ]) in panel (6) are the averaged values in 
the field of view, respectively. The dispersion of the optical thickness r = / a VQ ds and 
the lack of proportionality to the column density arise from the large dispersion of the 
absorption coefficient a u , which is calculated by the nonLTE level population (eq. [A5]) in 
the inhomogeneous molecular torus (see text). 
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Fig. 7. — Distribution of the integrated intensity of HCN(1 — 0) line overplotted by the 
contour of r > 1. The data are of y = 2 x 10 -9 . One can observe discordance of the bright 
regions (color scale) and the optically thick regions (see text). 



m 



0.5 

i 0.0 

-0.5 
-1.0 

















10'" 




1.0 






0.5 


10 w 


r-i 


0.0 








10" s 




-0.5 


10" 10 




-1.0 




lo ±l j o.u 




-15 -10 -5 f) 5 fO 15 
z|pc] 



■15 -10 -5 
z|pc| 



10 15 



-15 - ID -5 5 10 15 
z|pc| 



Fig. 8. — The distributions of Atq (lines with crosses) and the density (lines with diamonds) 
along three lines of sight indicated in Fig. 7. The lowest density («h 2 = 10~ 8 cm~ 3 ) 
seen in these panels is the cut-off density adopted in the hydrodynamic simulation. These 
three panels represent different circumstances: in panel (a) emission from a dense clump 
at z ~ 3 is seen through a tenuous ambient, in panel (b) the optical thickness is large due 
to accumulation of Ar of low density atmosphere encompassing the large scale height, and 
in panel (c) the intensity is strong because of the stimulated emission by the population 
inversion within a dense region (bh 2 < ftcrit) around z — 0, respectively. 
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Fig. 9. — The distribution of line profiles at the points with 10 pc spacing. Abscissa denotes 
the velocity component along the line of sight measured from the position of the torus. The 
binning size of V r is taken to be 16 km s _1 . In each panel, black lines are HCO + (l — 0) and 
red lines are HCN(1 — 0) lines. Multi-peak profiles at the points labeled (A) and (B) reflect 
the inhomogeneous structure of the torus. 
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Fig. 10. — The averaged line profiles over the whole field of view. The solid line denotes 
HCO + , and the dashed line denotes HCN profile. The inhomogeneous distribution of line 
profiles seen in Fig. 9 is compiled to a single Gaussian-like profile with width Av ~ 200 km 
s- 1 . 
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Fig. 11. — Density dependence of the ratio n\jn§ (the excitation temperature T ex ). Blue 
lines represent HCN, and black lines do HCO + transitions for various kinetic temperatures: 
long-dashed lines are T = 10K, three-dot-dashed lines are T = 30K, dot-dashed lines are 
T = 100K, short-dashed lines are T = 200K, and solid lines are T = 1000K. In both of 
the extremities of low and high density, ni/no for HCN and HCO + molecules coincide. See 
Appendix for details. 
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Fig. 12. — The fractional level population (fj) distributions of HCN rotational line for 
various densities in the optically thin limit with T = 100 K. In this panel, the dotted lines 
denote nn 2 = 10 5 cm -3 , the dashed lines tih 2 = 10 7 cm -3 , the dot-dashed lines nn 2 = 10 9 
cm -3 , and the solid lines the Boltzmann distribution. 
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Fig. 13. — The degree of population inversion 5 pop as a function of temperature and density 
for HCN(1 —0) (left panel) and HCO + (l — 0) (right panel) transitions which is calculated in 
the optically thin limit. Contours are from S pop = 1.0 to 1.8 with 0.2 spacing, and population 
inversion occurs within the outermost contour (gray regions). 



